ICM 2002 • Vol. Ill • 1-3 
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Abstract 

Multiple scale homogenization problems are reduced to single scale prob- 
lems in higher dimension. It is shown that sparse tensor product Finite Ele- 
ment Methods (FEM) allow the numerical solution in complexity independent 
of the dimension and of the length scale. Problems with stochastic input data 
are reformulated as high dimensional deterministic problems for the statistical 
moments of the random solution. Sparse tensor product FEM give a deter- 
ministic solution algorithm of log-linear complexity for statistical moments. 
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1. Introduction 

The numerical solution of elliptic problems with multiple scales in a bounded 
domain D can be achieved either by analytic homogenization [2], through 
asymptotic analysis or by specially designed Finite Element spaces to capture the 
fine scales of the problem 0, 0- In asymptotic theory of homogenization, the 
fine scale of the solution is averaged and lost in the homogenized limit. Fine scale 
information can be recovered from so called correctors which must be calculated 
separately. An alternative which we pursue here is to "unfold" the homogenization 
problem into a single scale problem in high dimension - so that homogenized and 
fine scale behaviour are still coupled 0. As we will show here, this limiting, high- 
dimensional single scale problem can be solved numerically by sparse tensor product 
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FEM in complexity comparable to the optimal one for single scale problems in the 
physical domain D. 

We also show how the same idea can be applied to the fast deterministic 
calculation of two and M point spatial correlation functions of random solutions to 
elliptic problems. 



2. Homogenization problem 

In a bounded domain D C R d with Lipschitz boundary T = dD 7 we consider 
the elliptic problem in divergence form 

-div(A e O)Vw e ) = f(x) in D, u e = on 3D . (2.1) 

Problem (|2.1|l has multiple separated scales in the sense that 

i £ (i)=i(i,{^J,ieJ} (2.2) 

where Y = (0, l) d denotes the unit cell and we denote by the unique element 
in 7L d such that x £ e([|] y + Y) and set {§ } Y ~ § - [f]y £ F. In J22J), the 
function A(x,y) £ L°°(D x is F-periodic with respect to y and satisfies, for 

every (x, y) £ D xY, and some < a < 1: 

G R d : a ICI 2 < e A(x, y) £ < a^] 2 . (2.3) 

Then problem Ij2.1|l admits, for every / £ L 2 (D), a unique weak solution: 

u s eH^(D): B"(u^v) = (f,v) V« £ H^D) (2.4) 

where B e (u,v) — J D Vw • A(x, ^)Vudx. As it is well known, as e — > 0, u e — > u° 
in L 2 (D) strongly and in H 1 (D) weakly, and m° is the solution of the homogenized 
problem 

-div(A°{x)Vu°) = f in D, u° = on dD , (2.5) 

and formulas for A°(x) are available jUEJ. The lack of strong iJ 1 (D)-convergence 
indicates that in the limit e — > 0, information on the fine scale of u £ is lost. It can 
be recovered by calculating so-called correctors by either differentiating u e or by 
solving a second problem of the type Q2.1JI . Both approaches are not attractive as 
basis for numerical solution methods: differentiating a numerical approximation of 
vP will reduce convergence rates, and solving l|2.1() for correctors amounts to solving 
a problem akin to the original one. 

A variant of the two-scale convergence, originally due to JOIQ, and recently 
developed in PJ, allows to obtain a single scale, "unfolded", limit problem which 
gives u° and essential information on oscillations of u s . To describe it, define for 
every ip £ L 2 (D) the "unfolding" operator 

\ / r T* "i \ 

:= { ' X ' ' > ^ (2.6) 
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for x e D,y 6 Y. Then there holds 0. 



Proposition 2.1. Assume that 



A e (x,y) :=%(A e )(x,y) -» A{x,y) a.e.{x,y) E D x Y. 



(2.7) 



Then there exists in u G H { 
satisfy 



(D) such that, as e — > ; the solutions u £ of \2.$ 



m(D) 



u Q eHl{D) 



(2.8) 



and there exists (j>(x,y) G L 2 (D , Hp Cr (Y) /Wl) such that, as e — > 0, 
%{V x u e ) — ^ V x u° + \7 y cf) in L 2 (L> x Y) d . 



L 2 (DxY) d 



(2.9) 



The functions u° G Hq(D), cj> G L 2 (L>, 7?p er (F)/R) sotoe ffce "unfolded" limiting 
problem: find u° G H^(D), cj> G L 2 (D, H^ eI (Y)/U) such that 



Remark 2.2. i) Problem (|2.1()(l is independent of e and therefore a single scale 
problem. 

ii) As e -> 0, u £ (a;) -> in L 2 (£>) and V^u 5 -> V^m + V y 0(-, -) strongly in 
L 2 (i?,R' i ). Therefore, u° coincides with the solution of the homogenized problem 
l|2.5|l and cj> retains information on the oscillations of u £ as £ — > 0. Information 
on gradients of m £ is introduced at the price of higher dimension of the limiting 
problem: while l|2.5|l is posed onDcR i , (|2.10() must be solved on D x Y C R 2d . 

iii) As H2.10f> is derived by weak convergence methods, no information about the 
boundary behaviour of u £ is obtained by solving Ij2.1t)|l . 

iv) In l|2.4l) . / = f(x) was assumed independent of e. Loadings of the type f(x, xje) 
may equally be accommodated. 

v) Property (|2.7H holds under certain conditions on A e , for example if A(x, y) = A(y) 
or if A G i^C (£>)). 

The limit problem (|2.1U|) is well-posed: 
Proposition 2.3. Assume Then there is C > such that 



/ / {V x v + V y ip} T A(x,y){V x u° + V y (f>}dydx 




D Y 



(2.10) 




D 



Vu G Hjj(D), G L 2 (D,H^ eI (Y)/U) : 
B(u,(t>;u,4>) > C(\\u\\ 2 H1(D) + \\<f>\\l 2(D!H ^ AY)) ) . 



(2.11) 
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The proof is immediate if we observe l|2.3[) and 
| V x u + V y 0| 2 = IV^uI 2 + 2\7 x u ■ 
and that, due to <j> G L 2 (D, if* (y)/]R), it holds 



V ' x u ■ V ' y <p dy dx — J V x u- J \7 y <pdy dx = . 
D Y D Y 

The limit problem (|2.10Jl admits the following regularity: if A(x, y), dD are smooth, 
then / G H~ 1+k (D) implies 

u° G H 1+k {D), <$> G H k (D, C^r(F)) • (2-12) 

3. Sparse finite element discretization 

We discretize the unfolded limiting problem H2.10|l by a sparse Finite Element 
Method (FEM) (e.g. [H] and the references there). To this end, assume D is a 
bounded Lipschitz polyhedron with straight faces and let {T/ 3 }, {T^} be sequences 
of nested, quasiuniform meshes in D resp. Y consisting of shape-regular simpliccs 
T of meshwidth hi = 0(2~ e ), and such that the periodic extension of beyond 
Y is regular. 

Let further p > 1 be a polynomial degree. Then we denote 

V"/ = Sq' 1 (D, Tf) = {u G Hl{D) : VT G T t D : u\ T G P P (T)} , 

V$ = S&{Y,T?) = {ue Hp Cr (Y) : VP G T/ : u\ T G P P (T)}/R, 

^,0 = SP-lfl( D)T D) = {ue L 2 [D) . yr £ T D . £ ^ p _ l(T)} 7 

where V P (T) denotes the polynomials of total degree at most p on T. 

Since the triangulations are nested, the Finite Element spaces are hierarchical: 

c r,y c-.-c v* 3 " c-.-c 

^ C V£ C-.-C V£ C-.-C Hl eT {Y). 

Define orthogonal projections P^ J : H^(D) -> V^ j , j = 0, 1 and P£ : H* ex (Y) 
7 Y and set P~ lj ' := 0, Pf 1 



and set P D lj := 0, P Y 1 := 0. Then we have the increment or detail-spaces 



wt? := i/-;,' /';/':i;r. *= 0,1,2,... 

:= (P^-P^ 1 )^, £ = 0,1,2,... 
and, for every L > 0, the multilevel decompositions 

Vy= j =0,1 (3.3) 

o<e<L o<e<L 
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which are orthogonal in fP(£>) resp. Hp CI (Y). 

For L > 0, define the full tensor product space 

S L = v£>°®Vf = W e j?®W^L 2 {D,Hl eT {Y)/U). (3.4) 

0<iJ'<L 

Assume the regularity (|2.12(1 with some k > 0. 

Then l|2.11|l shows that the finite element approximation 

{u L , 4> L ) G Vp 1 x S L : B(u L , 4> L ; v L , ^ L ) = (/, v L ) V(v L ^ L ) G y£' x x 5 L (3.5) 

exists, is unique and a tensor product argument show that it satisfies the asymptotic 
error bounds 



\u° - u L \\ H i {D) + ||0 - <t> L \\L*(D,m(Y)) 

< Ch™ a(p,k ' {\\u°\\ H k + i (D) + \\(f>\\ L i(D;H k + l (Y)) + W^WH^iD-H^Y))} 



We note that dim(l^" 1 ) = 0{2 Ld ) = dimV^, whereas dim(S L ) = 0(2 2Ld ) = 
0(dim(V^' 1 ) 2 ) as L — > oo, since S L is a FE-space in D x Y C R 2d . Also, note that 
in (|3.6|) the full regularity l|2.12|l of <f> was not used. 

The large number of degrees of freedom in S L due to high dimension renders 
the unfolded problem 12.10fl impractical for efficient numerical solution. A remedy 
is to use the sparse tensor product spaces and the references there) 

S L := W l jf ® W Y ' (3.7) 
o<e+e'<L 

which have substantially smaller dimension than the full tensor product spaces l|3.4|) . 

The joint regularity l|2.12|l in x and y of (f>(x, y) allows to retain the error 
bounds (|3~B|l with S L in place of S L : 

Proposition 3.1. Let (u L ,cj) L ) G V^' 1 x S L denote the sparse FE-solutions which 
satisfy \S. .5)) with S L in place of S L . 

Ifu° G H 1+k (D), <j) G H^D.H^+^Y)), for some k>l, then 



(3.8) 



- u L \\ H i {D) + ||0 - <P l \\l 2 (d,h^(y)) 

< C(L + l)i hf aM {\\u°\\ H ^ {D) + HWhhd-m^HY))} 
and the total number of degrees of freedom is bounded by 

dim(V^" 1 ) + dim(5 L ) < CL2 dL . (3.9) 

Remark 3.2. i) Since dim(V^' 1 ) < C2 dL , computation of (w,0 L ) in D x Y re- 
quires, up to L = 0(| logft-ij), the same number of degrees of freedom as the FE 
approximation of u° from (|2.5|l in D. Moreover, the FE approximation of l|2.10|l 
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does not require the determination of A (x) . In addition, the convergence rate (|3.8|l 
is, up to (L + 1)2, equal to the rate (|3.ti[) for the full tensor product-spaces. 

ii) Sparse tensor products of finite elements in one dimension were proposed by 
Zenger and his students in the 1990ies for the efficient solution of partial differential 
equations in three dimensions (see [5] and the references there). While allowing 
similar convergence rates as the full tensor product spaces, extra regularity of the 
solution (generally not available in non-smooth domains) is required. In . 7|> , sparse 
tensor products of standard finite element spaces in D resp. Y C R d are taken and 
realistic regularity of the solution available from the structure of the limit problem 
H2.I0f> was used. 

iii) The approach generalizes to problems with M > 2 scales. The unfolded problem 
is then posed on a product domain 

D xY 1 x ■■■ x Y M -i C U Md 

and an approximation of order (L+l)~ ^~ / l ™ ln (P' fc ) can b e obtained with 0(L M ~ X 2 dL ) 
degrees of freedom. 

iv) Explicit construction of the sparse tensor product-space S L requires bases for the 
detail- spaces W^ , W^. These are available, for example, via suitable semiorthog- 
onal wavelet bases (e.g. 0]). These bases allow also for optimal preconditioning of 
the linear system corresponding to l|3.5l) . 

v) We discussed here only the diffusion problem l|2.I|l . The results can be generalized 
to Elasticity, and the Stokes Equations. 

vi) If only partial periodicity or patchwise periodic patterns are present, the unfold- 
ing approach works equally. See |3] for details. 

4. Stochastic data 

Let (f2,£,P) be a cx-finite probability space and A(x) £ L°°(D, IR S3 ^q) satisfy 
for every £ G ]R d 

3a,P>Q:a\tf A(x)Z<l3\t\ 2 ■ (4-1) 
For a random source term / £ L 2 (fl, dP; L 2 (D)), consider the Dirichlet problem 

L(x, d x )u = -V • A(x) Wu(x, to) = f(x, w) in D, u = on dD (4.2) 

for P-a.c. uj £ CI. The random solution u(x, oS) of l|4.2|l is searched in the Bochncr- 
space 

H 1 {D):=L 2 {Q 1 dP;H^(D))=L 2 (n i dP)(g > H^(D). (4.3) 
We note that L 2 (Cl, dP), equipped with inner product 

(u,v) = J u(u)v(u)dP(ui), (4.4) 
n 
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is a Hilbert space. The variational form of (|4.2|l reads: find u 6 Hq(D) such that 
A(u,v) = (f,v) C 2 {D) VveHl(D) (4.5) 

where 

A(u, v) = ((A <g) id)(V ® id)u, (V ® id)v)&{Dy 

and (fiV)cz(D) — Je> (f(x, ■), -))<iz;. The form „4(-, •) in Ij4.5[l is coercive on 
Hl(D) x Hl(D), implying the existence of a unique random solution u(x,w) of 
(|4.5|l . Numerical solution of (|4.5|l that involves a FEM in D and Monte-Carlo in 
Q is prohibitively expensive. Alternatively, we might try to compute directly the 
statistics of u(x,uj). For example, the mean field E u (x) = J Q u(x,uj)dP(uj) solves 

L(x, d x )E u = E f in D, E u = on dD . (4.6) 
For x,x' € D, define the two point correlation function 

C u {x,x') = {u{x,-),u{x',-)). (4.7) 
Then the variance of u(x, lu) is given by 

(Var M (x, -)) 2 := (E u (x)) 2 - {C u {x,x)f . (4.8) 

The two-point correlation C u (x, x') is the solution of a deterministic problem in 
D x D. Formally 

L{x, d x ) L(x', d x ,)C u = C f in D x D , (4.9) 
and in variational form: 

C u S H$(D x D) : Q(C U , C v ) = {C f , C v ) VC* 6 H$(D x D) (4.10) 

where H^(D x D) := Hq(D) ® H$(D) and 

Q(C U ,C V )= J V xy C v A(x)®A(y)V xy C u dxdy,V X y - V ' x §\ . 

DxD 

Proposition 4.1. TVie two point correlation C u of the random solution u{x,uj) is 
the unique solution of the deterministic problem \4-lU\j in D x D . 

Hence to get two point correlation functions, Monte-Carlo can be traded for a 
deterministic problem in high dimensions. The key to its efficient numerical solution 
lies in the regularity of the solution: if the mean field problem H4.fif) satisfies a shift- 
theorem at order s > 0, i.e. E f G H- 1+S (D) =>• E u e H 1+S (D), then 

C f E H- 1+S '- 1+S (D xD)=^C u e Hl+ S ' 1+S {D x D) . (4.11) 
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A similar regularity result in weighted spaces of mixed highest derivatives holds if 
DcU 2 has corners ^Uj- This regularity in spaces of mixed highest derivatives 
allows to approximate C u from the sparse tensor product FE spaces 

V L := Wo^g/Wp 1 (4.12) 

0<t+i'<L 

of dimension dim(V i ) < CL2 dL at a near optimal convergence rate: if C£ G V L 
denotes the FE approximation of C u , it holds 

\\Cu - C^\\ m . 1(DxD) <C\logh L \i h? a(s ' p) \\C u \\ H , +1 ,. +HDxD) . (4.13) 

Moreover, using a semiorthogonal wavelet basis of the detail-spaces Wfj 1 in 13.21) . 
we can design an algorithm which computes to the order of the discretization 
error H4.13|l in 0(NlL a<1+2 ) operations where Nl — 0(2 dL ) denotes the number of 
degrees of freedom in D (see f° r details). 

As for homogenization problems with multiple scales, M point correlation 
functions of u(x,lu) can be approximated at the rate | log/ii|^ M_1 '/ 2 / l 1 ° m ^' p ^ ) with 
0(Nl L M ~ l ) degrees of freedom . Let us also remark that high regularity in l|4.13|l 
corresponds to strong spatial correlation of u(x, uj). 

In the spatially uncorrelated limit, formally C/(x, y) — S(x — y) and we have 
for smooth dD, A(x), d < 3: 

C u g H 1+S ' X+S (D x D), < s < 1 - j , 

so that only the low convergence rates | log/iij^ h L 4 e for follow. This is due 
to the singular support of C u being the diagonal {{x,y) : x = y}. The efficient 
approximation of such C u is the topic of ongoing research. 
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